from sympy import *
from sympy.matrices import Matrix
from sympy.physics.vector import ReferenceFrame, curl, divergence, gradient
import plotly.graph_objects as go
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
Curl and grad helper functions.
This code still needs testing functions written for it.
def nabla(tuplefield, kind='curl'):
'''
tuplefield: tuple
tuple of expressions involving x,y,z
'''
R = ReferenceFrame('R')
try:
x,y,z = sorted(list(tuplefield.free_symbols), key=lambda i: i.name)
except:
return NameError('x,y,z not defined')
# turn tuple into vector field using reference frame
field = rhopf.subs({x:R[0],y:R[1],z:R[2]})
field = field[0] * R.x + field[1] * R.y + field[2] * R.z
if kind == 'curl':
field = curl(field, R)
if kind == 'divergence':
field = divergence(field, R)
# turn vector field back into vector
field = field.subs({R[0]:x,R[1]:y,R[2]:z})
field = (field.dot(R.x), field.dot(R.y), field.dot(R.z))
return field
This is a class that performs a stereographic projection between $\mathbb{R}^n$ and $S^n$.
class StereographicProjection:
def __init__(self, special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3):
'''
special_point: tuple, default: (0,0,1)
special reference point on Sn from which to perform stereographic projection
sphere_centre: tuple, default: (0,0,0)
centre of the sphere in R^{n+1}
sphere_radius: float, default: 1
radius of the sphere
n: int, default: 2
dimension of the plane of stereographic projection
'''
self.sp = Matrix(special_point)
self.sc = Matrix(sphere_centre)
self.r = sphere_radius
self.n = n
# assert special point lies on sphere
assert (self.sp - self.sc).dot(self.sp - self.sc) == self.r
# assert special point does not lie on Rn
assert self.sp[-1] != 0
def project_Rn_to_Sn(self, v):
'''
vector on Rn expressed as an (n+1)-tuple with last entry null
'''
v = Matrix(v)
# assert point lies on Rn
assert len(v) == self.n + 1
assert v[-1] == 0
t = Symbol('t')
w = self.sp + t * (v - self.sp) - self.sc
t = solve(w.dot(w) - self.r, t)
if not t: return nan
return tuple(simplify(self.sp + t[1] * (v - self.sp)))
def project_Sn_to_Rn(self, v):
'''
vector on Sn expressed as an (n+1)-tuple
'''
v = Matrix(v)
# assert point lies on Sn
test = N(v.dot(v))
if not test.free_symbols:
assert np.isclose(float(test), 1)
t = self.sp[-1] / (self.sp[-1] - v[-1])
v = self.sp + t * (v - self.sp)
if not test.free_symbols:
return tuple([float(i) for i in N(v)])
return tuple(simplify(v))
This is the version of the Hopf map that is formulated using the Riemann sphere.
The Riemann_Hopf_Map function returns a vector field by performing the following series of maps:
$\mathbb{R}^3 \xrightarrow{\text{stereo}} S^3 \simeq \mathbb{C}^2 \xrightarrow{z_1/z_2} \mathbb{C} \simeq \mathbb{R}^2 \xrightarrow{\text{stereo}} S^2$.
def Riemann_Hopf_Map():
# consider any point (x,y,z) in R3
x,y,z = symbols('x y z', real=True)
# stereographic projection to S3
stereo1 = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
a,b,c,d = stereo1.project_Rn_to_Sn((x,y,z,0))
# compute the Riemann Hopf map by dividing two complex numbers
w = (a + I*b) / (c + I*d)
rew = simplify(re(w))
imw = simplify(im(w))
# stereographic projection to S2
stereo2 = StereographicProjection(special_point=(0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
return factor(stereo2.project_Rn_to_Sn((rew,imw,0)))
rhopf = Riemann_Hopf_Map(); rhopf
The preimage of a point under the Hopf map is a closed loop. The preimages corresponding to any two distinct points are a pair of linked closed loops.
We can explicitly compute the preimages of the Hopf map by noticing $$\eta(z_1,z_2)=\dfrac{z_1}{z_2} = \dfrac{r_1}{r_2}e^{i(\theta_1-\theta_2)} = R e^{i\phi}$$ And since $r_1^2 + r_2^2 = 1$, we have that $$z_1 = \sqrt{\dfrac{1}{R^2+1}}e^{i(\phi+t)},\quad z_2 = \sqrt{\dfrac{R^2}{R^2+1}}e^{it},$$ where $t$ is an arbitrary real parameter. Clearly $(z_1, z_2)$ is $2\pi$ periodic in $t$ and therefore forms a closed loop.
def Inverse_Riemann_Hopf_Map(vector):
'''
Given any tuple (x,y,z) in S2, computes the preimage in R3 under the Riemann Hopf Map.
'''
# map S2 to R2 equivalent to C
stereo1 = StereographicProjection(special_point=(0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
a, b, _ = stereo1.project_Sn_to_Rn(vector)
# compute preimages of Hopf map in C2
# Rewrite a + bi in polar form R * exp(i * Phi)
R = sqrt(a**2 + b**2)
if N(a).free_symbols or N(b).free_symbols:
Phi = atan(b/a)
elif np.isclose(float(abs(b / a)), float(pi / 2)):
Phi = [-1, 1][b / a >= 0] * pi / 2
else:
Phi = atan(b/a)
# Compute radii of preimages
r_1 = R / (R**2 + 1) ** .5
r_2 = 1 / (R**2 + 1) ** .5
# Write parametrised preimages
t = Symbol('t')
vector = (r_1 * cos(Phi + t), r_1 * sin(Phi + t), r_2 * cos(t), r_2 * sin(t))
# map C2 to R3
stereo2 = StereographicProjection(special_point=(0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
a, b, c, _ = stereo2.project_Sn_to_Rn(vector)
return (a, b, c)
def compute_hopf_preimages(vector):
'''
Compute a set of Hopf map preimage points for a tuple in R3.
'''
vector = Inverse_Riemann_Hopf_Map(vector)
t, theta = symbols('t theta')
vectorx = lambdify((t, theta), vector[0], 'numpy')
vectory = lambdify((t, theta), vector[1], 'numpy')
vectorz = lambdify((t, theta), vector[2], 'numpy')
mesht = np.mgrid[0:2*np.pi:300j]
meshtheta = np.mgrid[0:2*np.pi:300j]
x = vectorx(mesht, meshtheta)
y = vectory(mesht, meshtheta)
z = vectorz(mesht, meshtheta)
return (x,y,z)
x, y, z = compute_hopf_preimages((1,0,0))
a, b, c = compute_hopf_preimages((sqrt(2)/2, 0, sqrt(2)/2))
fig = go.Figure()
fig.add_trace(go.Scatter3d(
x=x, y=y, z=z,
marker=dict(
size=0,
),
line=dict(
color='darkblue',
width=2
)
))
fig.add_trace(go.Scatter3d(
x=a, y=b, z=c,
marker=dict(
size=0,
),
line=dict(
color='darkred',
width=2
)
))
fig.update_layout(
template='ggplot2',
scene = dict(
xaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
yaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
zaxis = dict(
showbackground=False,
showticklabels=False,
),
),
width=700,
margin=dict(
r=10, l=10,
b=10, t=10)
)
fig.show()
The linking number of the preimages of any two distinct points for the standard Hopf map is a conserved quantity. This is formalised through the the Hopf invariant which can be defined as $$H(f) = \int_{S^3}\omega\wedge d\omega$$ where $\alpha$ is a generator of $H_{\text{DR}}^2(S^3)$ and $\omega$ is a $1$-form such that $f^*\alpha=d\omega$ and $f$ is a map from $S^3$ to $S^2$.
Programmatic implementation of this integral coming soon...
The Hopf texture can be visualised using the PT construction. We set a plane of orientation, say $n_z=0$, and plot the surface corresponding to vectors pointing along this orientation coloured by their phase.
def compute_hopf_texture():
'''
Draws the PT construction for the Riemann Hopf Map.
'''
# declare coordinate symbols and parametric symbols
x,y,z = symbols('x y z')
t, theta = symbols('t theta')
# compute the preimages of the z = 0 orientation
vector = Matrix(Inverse_Riemann_Hopf_Map((x,y,0))).subs({x**2+y**2:1, atan(y/x):theta})
vectorx = lambdify((t, theta), vector[0], 'numpy')
vectory = lambdify((t, theta), vector[1], 'numpy')
vectorz = lambdify((t, theta), vector[2], 'numpy')
mesht, meshtheta = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j]
x = vectorx(mesht, meshtheta)
y = vectory(mesht, meshtheta)
z = vectorz(mesht, meshtheta)
# alongside coordinates return the orientation modulo pi as a colour scale
return x, y, z, 2 * (meshtheta % np.pi)
x, y, z, colour = compute_hopf_texture()
fig = go.Figure(
go.Surface(x=x, y=y, z=z, surfacecolor=colour, colorscale='greys'),
)
fig.update_layout(
template='ggplot2',
scene = dict(
xaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
yaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
zaxis = dict(
showbackground=False,
showticklabels=False,
),
),
width=700,
margin=dict(
r=10, l=10,
b=10, t=10)
)
fig.show()
The corresponding twist, $\mathbf{n} \cdot \nabla \times \mathbf{n}$, can be visualised through its isosurfaces.
twist = factor(Matrix(nabla(Matrix(rhopf), kind='curl')).dot(Matrix(rhopf))); twist
x,y,z = symbols('x y z')
numpy_twist = lambdify((x,y,z), twist, 'numpy')
X, Y, Z = np.mgrid[-5:5:40j, -5:5:40j, -5:5:40j]
values = numpy_twist(X,Y,Z)
fig = go.Figure(
data=go.Isosurface(
x=X.flatten(),
y=Y.flatten(),
z=Z.flatten(),
value=values.flatten(),
opacity=0.9,
surface_fill=0.4,
isomin=-1,
isomax=1,
surface_count=3,
caps=dict(x_show=False, y_show=False),
colorscale='Cividis',
),
)
fig.update_layout(
scene = dict(
xaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
yaxis = dict(
showbackground=False,
showticklabels=False,
showaxeslabels=False,
),
zaxis = dict(
showbackground=False,
showticklabels=False,
),
),
width=700,
margin=dict(
r=10, l=10,
b=10, t=10)
)
fig.show()
The pitch of the vector field describes the distance and axis over which a full rotation occurs. It can be defined as the positive eigenvector of $$\Pi_{ij} = \dfrac{1}{4}\epsilon_{ilk}[n_l \partial_k n_j + n_l \partial_j n_k - n_j n_l n_m \partial_m n_k] + \dfrac{1}{4}\epsilon_{jlk}[n_l\partial_k n_i + n_l \partial_i n_k - n_i n_l n_m \partial_m n_k]$$
Still in progress
The vector field can be explicitly visualised. It is uniform in the far field with a knotted localised structure near the origin.
x, y, z = symbols('x y z')
# compute the preimages of the z = 0 orientation
vectorx = lambdify((x, y, z), rhopf[0], 'numpy')
vectory = lambdify((x, y, z), rhopf[1], 'numpy')
vectorz = lambdify((x, y, z), rhopf[2], 'numpy')
x, y, z = np.mgrid[-5:5:10j, -5:5:10j, -5:5:10j]
u = vectorx(x, y, z)
v = vectory(x, y, z)
w = vectorz(x, y, z)
fig = plt.figure(figsize=(5,5))
ax = fig.gca(projection='3d')
ax.quiver(x, y, z, u, v, w, length=1)
plt.show()
We can compare this unusual twist structure to the free energy. In the single constant approximation, the Frank free energy density reduces to $$\dfrac{1}{2}K(\nabla \mathbf{n})^2 + Kq_0\mathbf{n}\cdot(\nabla \times \mathbf{n}),$$ where the constant can now be neglected.
This section is coming soon once I know what I'm doing.
This is the version of the Hopf map formulated using quaternions. It connects the map to the rotation group.
The Quaternion_Hopf_Map function returns a vector field by performing the following series of maps:
$\mathbb{R}^3 \xrightarrow{\text{stereo}} S^3 \xrightarrow{qpq^*} \mathbb{H} \simeq S^2$.
In the interest of time, I might abandon this section.
def Quaternion_Hopf_Map(x,y,z):
pass
These are test functions that should all pass if the code is working properly.
def test_project_Rn_to_Sn():
# R2 to S2
stereo = StereographicProjection(special_point = (0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
# general point
vector = (4,5,0)
output = stereo.project_Rn_to_Sn(vector)
assert output[0] / (1 - output[-1]) == vector[0]
# origin
vector = (0,0,0)
output = stereo.project_Rn_to_Sn(vector)
assert output == (0,0,-1)
# circle point
vector = (0,1,0)
output = stereo.project_Rn_to_Sn(vector)
assert output == vector
# R3 to S3
stereo = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
# general point
vector = (4,5,1,0)
output = stereo.project_Rn_to_Sn(vector)
assert output[0] / (1 - output[-1]) == vector[0]
# origin
vector = (0,0,0,0)
output = stereo.project_Rn_to_Sn(vector)
assert output == (0,0,0,-1)
# circle point
vector = (0,1,0,0)
output = stereo.project_Rn_to_Sn(vector)
assert output == vector
# symbols
x,y,z = symbols('x y z')
vector = (x,y,z,0)
output = stereo.project_Rn_to_Sn(vector)
assert output[-1] == (x**2 + y**2 + z**2 - 1)/(x**2 + y**2 + z**2 + 1)
test_project_Rn_to_Sn()
def test_project_Sn_to_Rn():
# R2 to S2
stereo = StereographicProjection(special_point = (0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
# general point
vector = (0.5, 0.5, sqrt(2) / 2)
output = stereo.project_Sn_to_Rn(vector)
assert np.isclose(output[-1], 0)
assert np.isclose(output[0] / output[1], float(vector[0] / vector[1]))
assert np.isclose(1 - float(vector[0]) / output[0], float(vector[-1]))
# origin
vector = (1,0,0)
output = stereo.project_Sn_to_Rn(vector)
assert output == vector
# symbols
x, y = symbols('x y')
vector = (x, y, x * y)
output = stereo.project_Sn_to_Rn(vector)
assert output[-1] == 0
assert output[0] / output[1] == vector[0] / vector[1]
assert 1 - vector[0] / output[0] == vector[-1]
# R3 to S3
stereo = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
# general point
vector = (0.5, 0.5, 0, sqrt(2) / 2)
output = stereo.project_Sn_to_Rn(vector)
assert np.isclose(output[-1], 0)
assert np.isclose(output[0] / output[1], float(vector[0] / vector[1]))
assert np.isclose(1 - float(vector[0]) / output[0], float(vector[-1]))
# origin
vector = (1,0,0,0)
output = stereo.project_Sn_to_Rn(vector)
assert output == vector
# symbols
x, y, z = symbols('x y z')
vector = (x, y, z, x * y * z)
output = stereo.project_Sn_to_Rn(vector)
assert output[-1] == 0
assert output[0] / output[1] == vector[0] / vector[1]
assert 1 - vector[0] / output[0] == vector[-1]
test_project_Sn_to_Rn()